function twoshiftpddyneps(N,D1,D2,M,amp)
% twoshiftpddyneps(N,D1,D2,M,amp)
% 	Parameters:
% 		N:	number of data points to include in time series
% 		D1:	first number of seconds to displace time series over
% 		D2:	second number of seconds to displace time series over
% 		M:	kernel size
% 		amp:	bound for the absolute value for sequence of perturbations
% 	Outputs:
% 		Only the plots of this code produces.
% 	Purpose of this code:  This code is a modification of the earlier code twoshiftpdeps.  Here the perturbations at each time step are dynamic making
% 	making the delay a function of time, but also the perturbation of the delay.  The absolute value of these perterbations is bounded above by the 
% 	value amp.


hold off;
eps = zeros(1,N);
eps (2:N) = randeps(N-1,amp);
ts = whitenoise(N,'w',1);
%sig = testsine(1:N);
%ts2 = ts+sig;
ts = lowpass(ts);
%ts2 = lowpass(ts2);
[psorig,freq] = mypsd(ts);
%% shift timeseries by S1 and then by S2
F11 = FDtesteps(ts,D1,M,@sinc);
F12 = FDtesteps(ts,D1,M,@sincdir);
F13 = FDtesteps(ts,D1,M,@sinc2dir);


epstot = cumsum(eps);
for j = 1:N
    delta(j) = epstot(j);
    prop = 0;
    delta(j) = delta(j).*prop;
end


for j = 1:N
    tsint10(j) = F11(j);
    tsint11(j) = F11(j) + delta(j).*F12(j);
    tsint12(j) = F11(j) + delta(j).*F12(j) + (1./2).*(delta(j)^2).*F13(j);
end

for j = 1:N
    Dadj = D1 - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end

tsint10 = tsint10./n0;
tsint11 = tsint11./n1;
tsint12 = tsint12./n2;

tsd1 = tsint11;

F21 = FDtesteps(tsd1,D2,M,@sinc);
F22 = FDtesteps(tsd1,D2,M,@sincdir);
F23 = FDtesteps(tsd1,D2,M,@sinc2dir);

for j = 1:N
    delta(j) = epstot(j);
    delta(j) = delta(j).*(1-prop);
end


for j = 1:N
    tsint20(j) = F21(j);
    tsint21(j) = F21(j) + delta(j).*F22(j);
    tsint22(j) = F21(j) + delta(j).*F22(j) + (1./2).*(delta(j)^2).*F23(j);
end

for j = 1:N
    Dadj = D2 - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end

tsint20 = tsint20./n0;
tsint21 = tsint21./n1;
tsint22 = tsint22./n2;

tsd2 = tsint21;

ps = mypsd(tsd2);

D = D1+D2;

F1 = FDtesteps(ts,D,M,@sinc);
F2 = FDtesteps(ts,D,M,@sincdir);
F3 = FDtesteps(ts,D,M,@sinc2dir);

for j = 1:N
    delta(j) = epstot(j);
end

for j = 1:N
    tsint0(j) = F1(j);
    tsint1(j) = F1(j) + delta(j).*F2(j);
    tsint2(j) = F1(j) + delta(j).*F2(j) + (1./2).*(delta(j)^2).*F3(j);
    
end
for j = 1:N
    Dadj = D - (M-1)/2; 
    Dint = round(Dadj); 
    k = (-(M-1)/2):1:((M-1)/2); 
    s = (Dadj-Dint)-k; 
    w = getWindow('blackman3',M);
    p = sinc(s).*w;
    q = sincdir(s).*w;
    q2 = sinc2dir(s).*w;
    
    intf0 = (p);% + eps.*q + (1./2).*(eps).^2.*q2);
    intf1 = (p + delta(j).*q);% + (1./2).*(eps).^2.*q2);
    intf2 = (p + delta(j).*q + (1./2).*(delta(j)).^2.*q2);
    %intf = (p + q);
    n0(j) = sum(intf0);    
    n1(j) = sum(intf1);    
    n2(j) = sum(intf2);
end
tsint0 = tsint0./n0;
tsint1 = tsint1./n1;
tsint2 = tsint2./n2;

tsonce = tsint1;
psa = mypsd(tsonce);

%% First Plot
%loglog(1:length(psorig),psorig,'b');
%hold on;
%err = abs(ps - psa)./abs(psa);
%xnn = abs(ps - psa);
%supp = sqrt(psa./xnn);
%loglog(1:length(ps),err,'g');
%loglog(1:length(ps),supp,'r');
err2 = abs(tsd2-tsonce);
errps = mypsd(err2);
%loglog(freq,errps,'g');
%hold on;
suppr = sqrt(psorig./errps);
figure;
loglog(freq,suppr,'b');



